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ABSTRACT 

Computations of drops collisions, coalescence, and other problems involving drops are presented. The 
computations are made possible by a finite difference/front tracking technique that allows direct solutions of the 
Navier-Stokes equations for a multi-fluid system with complex, unsteady internal boundaries. This method has been 
used to examine the various collision modes for binary collisions of drops of equal size, mixing of two drops of 
unequal size, behavior of a suspension of drops in linear and parabolic shear flows, and the thermal migration of 
several drops. The key results from these simulations are reviewed. Extensions of the method to phase change 
problems and preliminary results for boiling are also shown. 


INTRODUCTION 

The presence of a free surface that is not constrained to be more or less flat due to the action of gravity is 
perhaps the most important aspect of fluid flow in microgravity. The absence of gravity generally makes surface 
tension effects important at much larger length scales than we are used to on earth. The large amplitude surface 
motion possible when gravity is small or absent is generally highly nonlinear and thus difficult to analyze by 
conventional means. Such surface motion is, nevertheless, to be expected in microgravity environment and it is 
necessary to understand it to be able to work with free surface systems. Experiments in microgravity are expensive, 
and usually difficult as well. It is therefore desirable to be able to predict the evolution of free surface systems 
numerically. Numerical simulations can replace experiments in some cases and complement experiments in other 
cases by providing information that are difficult to measure. Numerical simulations can also aid in the design of 
space experiments. 

Numerical simulations of free surface flows have remained one of the frontiers of computational fluid dynamics since 
the beginning of large scale computations of fluid flow. Progress has however been much slower than for 
homogeneous flows and numerical simulations have not played the same role in multifluid and multiphase research 
as they have done for turbulence research, for example. Recently, however, a number of investigators have made 
considerable progress. We have been conducting simulations of multifluid and multiphase systems for several years 
and this paper reviews briefly those aspects of our work that are of particular relevance to microgravity fluid systems 
We start by a description of the numerical methodology, since it is critical to the success of our work, and then 
review a few applications. 

NUMERICAL METHOD 

The numerical method is based on writing one set of equations for the entire computational domain, 
independently of how many different fluids are involved. This is possible by allowing for different material properties 
in the formulation and adding singular terms at the phase boundaries to ensure that the correct boundary conditions 
are satisfied. The resulting "one-field" Navier-Stokes equations are: 

+ V • puu ~ -Vp + / + V fi(yu + V« r ) + J F s 5[x- Xy'jda 

Here, u is the velocity vector, p the pressure, and p and p are the discontinuous density and viscosity fields, 
respectively. / is a body force that can be used to initiate the motion. The surface forces, F a act only on the 
interface between the different fluids and appears in the current formulation multiplied by a three-dimensional delta 
function, S . The integral is over the entire front. It is important to note that this equation contain no 
approximations beyond those in the usual Navier-Stokes equations. In particular, it contains implicitly the proper 
stress conditions for the fluid interface. The momentum equation is supplemented by an equation of mass 
conservation, which for incompressible flows is simply 

V’M = 0 
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Combining this equation with the momentum equation leads to an elliptic equation for the pressure. Since the 
density and the viscosity are different for the different fluids, it is necessary to track the evolution of these fields by 
solving the equations of state which simply specify that each fluid particle retains it original density and viscosity: 

^- + u-Vp = 0; Hp + u V/i= 0 
at dt 

The momentum equation is discretized on a regular staggered grid using second order, centered finite differences for 
the spatial derivatives and a second order time integration scheme. The continuity equation, when combined with the 
momentum equation results in a pressure equation that is not separable as for homogeneous flow and is solved by a 
multigrid package (MUDPACK from NCAR). To advect the material properties, and to evaluate the surface tension 
term in the momentum equation, we track the interface between the different phases explicitly by using a moving 
grid of lower dimension than what we use for the conservation equations. This grid is usually referred to as a front. 
The one-field formulation used here is common to other techniques for multifluid flows such as the VOF (Volume of 
Fluid) and level set methods. In these methods, however, the phase boundary is not tracked explicitly, but 
reconstructed from a marker function. Explicitly tracking the interface avoids the difficulty of advecting such marker 
function and allows accurate evaluation of surface forces. 

Computing surface tension accurately is one of the most difficult part of methods intended for simulations of 
multifluid flows. Our current algorithm, which appear to be very satisfactory, is based on computing directly the 
force on each element by 

F a = <j> on x ids 
Jelem 

in a semi-implicit way. Here t is a tangent to the boundary of the surface element, n is the surface normal, and a 
is the surface tension coefficient. By computing the surface forces this way, we explicitly enforce that the integral 
over any portion of the surface gives the right value, and for closed surfaces, in particular, we enforce that the 
integral of surface tension forces is zero. This is particularly important for long time simulations where a failure to 
enforce this constrain can lead to unphysical motion of bubbles and drops. 

Since the boundary between the fluids (the front) usually undergoes considerable deformation during each run, it is 
necessary to modify the surface mesh dynamically during the course of the computations. The surface mesh is an 
unstructured grid consisting of points that are linked by elements. Both the points and the elements are arranged in a 
linked list, so it is relatively easy to change the structure of the front, including adding and deleting points and 
elements. Topological changes, such as when bubbles coalesce or drops break in two can also be accomplished by 
minimum effort. This is usually considered a major difficulty in implementing methods that explicitly track the 
front, but we have shown that with the right data structure these tasks become relatively straight forward. Although 
topology changes are easily done from a programming point of view, the physics is far from trivial. In reality, drops 
bounce off each other if the time when the drops are close is shorter than the time it takes to drain the film. Usually 
the film becomes very thin before it breaks and it would require excessive grid refinement to resolve the draining 
fully. At the moment we are dealing with this issue in a rather ad hoc way by simply changing the topology of the 
front at a prescribed time. However, considerable analytical work has been done on film draining and rupturing and 
we are currently exploring the possibility of combining such a model with our simulations. 

The method has been implemented for two- and three-dimensional flows on regular grids and for axisymmetric 
geometries using stretched grids to allow local grid refinement. It has been applied to a number of multifluid 
problems and tested and validated in a number of ways, not only to check the implementation, but also to assess its 
accuracy. Those tests include comparisons with analytical solutions for simple problems, other numerical 
computations, and experiments. The actual resolution requirement varies with the parameters of the problem. High 
Reynolds numbers, for example, generally require finer resolution than lower ones, as in other numerical 
calculations. We have also found that for problems where the surface tension varies, such as for contaminated 
bubbles and drops moving by thermal migration we generally require finer resolution than for flows where the 
surface tension is constant. However, in all cases have we found that the methods converges rapidly under grid 
refinement, and in those cases where we have other solutions we have found excellent agreement, even for modest 
resolutions. The method has also been extended in various ways to allow for simulations of heat transfer, surfactant 
effects, thermal migration, and phase changes. These effects can also be dealt with using the one-field formulation, 
thus allowing for a very efficient solution process. Examples of these validations are contained in various papers and 
dissertations listed at the end of theis paper. Our work on other physical systems, such as bubbly flows are also 
discussed in these references. 


536 



RESULTS 


We have examined the collision and coalescence of two drops in considerable detail. The results of our 
investigation of the head-on collision of drops can be found in Nobari, Jan, and Tryggvason (1996). In this paper we 
examine how the collision modes changed as a function of Weber number, and how the exact time of coalescence, 
where the film between the drops ruptures, affected the outcome of the collision. Figure 1 shows one simulation 
from this paper. This work has now lead to a collaboration with Professor C.K. Law of Princeton University who 
has conducted extensive experimental studies of drop collisions and a detailed comparison between his experimental 
results and the simulations has shown excellent agreement. The interplay between experiments and simulations has 
also greatly improved our understanding of the boundaries between the various collision regimes. Three-dimensional 
simulations of off-axis collisions are described in Nobari and Tryggvason (1996). The results there show good 
agreement with experimental investigations, but the study is not yet as detailed as for the head-on collisions. For 
drops of unequal size, we have examined the coalescence of initially stationary drops that are driven together by 
surface tension. The results, which have been submitted for publication (Nobari and Tryggvason, 1994), show that 
considerable mixing can take place due to vorticity injected into the large drop by the smaller one. 

With the collision of drops well under control, we have examined more complex problems. These include both large 
scale simulations of many drops as well as problems where additional physical effects have to be accounted for, such 
asheat transfer, thermal migration of drops, and phase changes. For isothermal drops we have examined the behavior 
of many drops in wall driven and pressure driven shear flows. In normal gravity such situations occur frequently in 
emulsions, for example, where the density of the drops is comparable with the continuous fluid. In microgravity, 
there is no restriction on the density difference and for finite Reynolds number flows we can expect significant effect 
due to the inertia of the drops. This investigation also has relevance to high pressure sprays where there are strong 
interactions between the drops and the continuous phase. In figure 2 we show several frames from a fully three- 
dimensional simulation of twelve drops in a periodic channel. The domain simulated is resolved by a 64^ grid and 
the simulations are carried out for sufficiently long time so the drop distribution reaches an approximately steady 
state. Using both three-dimensional simulations with a models number of drops and two-dimensional simulations 
with larger number, the effect of the various governing parameters has been explored. The simulations generally 
show that the drops migrate away from walls, forming a denser core near the centerline. A single drop, however, 
generally takes up a position approximately half way between the wall and the centerline (the so-called Segre-Silberg 
effect) and for dilute distributions we find local maximum in the concentration at this position. 

When gravity is absent, other effects, such as variable surface tension can dominate the dynamics. While thermal 
migration, where surface tension variations are brought about by variable temperature field, has been extensively 
studied for a single bubble or a drop, little is known about the interactions of many drops, particularly at finite 
Reynolds and Marangoni numbers. Figure 3 show an example of a fully three-dimensional simulation of two drops 
on a 32 by 64 by 1 28 grid. The drop surface and the temperature field in the centerplane are shown. The top and 
bottom boundaries are rigid walls at constant temperature, but the horizontal boundaries are periodic. This 
simulation, along with several two-dimensional ones of several drops show that the drops tend to form layers 
perpendicular to the temperature gradient. The simulations also show that even though surface tension is relatively 
low, the drops generally do not deform to any significant degree as they move. 

The methodology used for the problems discussed above can be extended to much more complex physics in a 
relatively straight forward way. Application to the relatively simple problem of solidification of pure material in the 
absence of flow is discussed in Juric and Tryggvason (1995) where validation studies for simple problems and 
applications to the growth of dendrites can be found. Extension to binary alloys where it is necessary to follow the 
composition of the material can be found in Juric (1996). For full multiphase problems where fluid flow, heat 
transfer, and phase change are taking place simultaneously, we combine the method for the energy equation used for 
the solidification problems with the method for the Navier-Stokes equations. The resulting code has allowed us to 
simulate several boiling problems and in figure 4 we show a simulation of the growth of a vapor bubble in 
superheated fluid. The domain, which is resolved by a 100 by 100 grid, is open at the top, allowing fluid to escape 
as the bubble grows. The bubble shape is perturbed slightly at time zero and as the bubble expands, these 
perturbations grow. The bubble shape at time zero and several subsequent times is shown in the frame to the left. 
The right frame shows the bubble shape and temperature contours at an intermediate time. A simplified version of 
this method has also been used to simulate solidification of drops colliding with a wall in the case where the density 
of the melt and the solid are the same. 
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CONCLUSIONS 


We have discussed a numerical methodology to simulate multifluid and multiphase flows in microgravity 
and shown a few application of these methods. These simulations have already lead to a better understanding of some 
aspects of drop collisions and coalescence and also demonstrate the versatility of the method. Overall, it seems 
reasonable to state that the fluid problem is under good control with the exception of how to handle the rupture of 
thin films, although our current ad hoc strategy seems to work well in some cases. In most cases the fluid motion 
is only a part of an engineering application and heat transfer and phase change must be dealt with in order to 
establish the understanding needed for a completely predictive capability. We have taken preliminary steps in this 
direction and examined the thermal migration of drops and various phase change problems. The capability to examine 
the detailed evolution of boiling fluids, for example, promises to yield a dramatic new insight and understanding. 
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